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ABSTRACT 



The objective of this thesis is to develop segmentation methods for multichannel 
and single channel images, and compare these methods. The segmentation algorithms 
are based on a linear model for the image textures and on inverse filtering to estimate 
the image textures and their regions. Two specific methods are compered 1) A 
multichannel filtering algorithm that simultaneously models the three separate signals 
representing the intensity of red, green, and blue as a function of spatial position and 
,2) A single channel model applied to a combined image resulting from performing a 
Karhunen-Loeve transformation on the three signal components. Results of the 
multichannel image segmentation and the Karhunen-Lbeve transformed one-channel 
image segmentation are presented and compared. 
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1. INTRODUCTION 



Segmentation techniques are among the most important considerations in the 
development of the automated image processing systems. Two related algorithms using 
2-D linear prediction models and the Karhunen-Loeve transformation for multichannel 
and color image segmentation are developed and compared in this thesis. 

The purpose of segmentation is to partition an image into a set of simpler 
homogeneous regions. The regions may consist of different gray level, different 
textures, colors, etc. In some cases an " image " may consist of several spectral 
components. For example, a color image consists of three separate signals representing 
the intensity of red, green, and blue as a function of spatial position. If we represent 
each of these signals by functions (n,m), Fj^ (n,m), and Fg (n,m), the image is 
represented by a vector quantity 



F (n,m) = 



'F^(n,m) 

Fb 

.Fg(n,m) , 



( 1 . 1 ) 






where n and m represent spatial coordinates. We call such an image, consisting of 
more than one two dimensional (2-D) signal, a multichannel image. 

In this thesis, a method based upon linear prediction is evaluated experimentally. 
This method has been developed [Ref. 1] for monochrome images and extended to 
color images [Ref. 2]. That method uses maximum likelihood {ML) and maximum a 
posteriori {MAP) estimation to segment multichannel images into regions of similar 
textures. The linear prediction is a filtering of the multichannel image to estimate the 
gray level at a particular spatial coordinate from the gray levels at neighboring 
positions. It is implemented as a 2-D linear filtering operation. The algorithm uses a 
previously-determined set of parameters corresponding to the mean of the data in each 
channel, the covariance matrix of the prediction error, and the weighting coefficients of 
the estimation filter for each specific texture type. 

The method discussed above is compared to a variant of this method based on 
the Karhunen-Loeve {K-L) transformation. The K-L transformation allows the several 
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components of an image to be combined into a single image that retains most of the 
energy in the original image. Hunt and Kubler [Ref 3] found that for image 
restoration, Karhunen-Loeve transformation followed by single channel image 
processing worked nearly as well as multichannel image processing. It was desired to 
see if the Karhunen-Loeve transformation would be equally effective for segmentation. 
In this part of the work, the K-L transformation has been developed to reduce the 
3-channel color problem to a 1-channel problem and the segmentation was performed 
for the one-channel image. Karhunen-Loeve transformation is based upon the 
statistical characteristics of an image. The advantage of this approach is 
computational savings; only about one ninth the number of computations is required 
for this method. 

The remainder of this thesis is organized as follows. Chapter II discusses the 
model and the algorithms used to perform the multichannel image segmentation 
employing techniques of linear prediction [Ref 4]. A general class of linear filtering 
models for texture is first presented. An algorithm is then developed to estimate the 
filter parameters from a multichannel image. Then, the multichannel image 
segmentation algorithm is described. 

Chapter III presents the models and the algorithms to perform the Karhunen- 
Loeve transformation and one-channel image segmentation. First, the algorithms for 
determining the eigenvectors and the eigenvalues of the correlation matrix are 
developed. Then, the transformation using a 3-channel image is presented. Finally, 
the one-channel image segmentation algorithm is discussed. The examples 
demonstrating the application of the segmentation methods for color images are 
presented and compared in Chapter IV. Chapter V has the conclusions about the 
multichannel image segmentation and one-channel image segmentation. 

In Appendix A, the Relaxation Method is described briefly as an alternative to 
the maximum a posteriori region estimation for monochrome images. Results are 
compared with the MAP method. In each of Appendices B and C, the description and 
use of the computer programs for the multichannel image segmentation and the one- 
channel image segmentation are presented. The computer program for multichannel 
image segmentation is contained in Appendix D. The Karhunen-Loeve transformation 
and one-channel image segmentation computer programs are contained in Appendix E. 
These programs are written in FORTRAN, compiled using Version 4.2 under the VAX 
/ VMS Version 4.1 operation system. 
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II. IMAGE SEGMENTATION USING A MULTICHANNEL FILTERING 

MODEL 



In this chapter, a multichannel image segmentation algorithm based upon a 2-D 
linear filtering model is presented. The multichannel images used in this work are color 
images with three channels representing the red, green, and blue components. The 
linear filtering model is used to develop approximate expressions for the multivariate 
probability density functions in terms of the filter error residuals for the entire set of 
points representing an image. The density expressions are used in the formulation and 
solution of the multichannel image segmentation problem. It is assumed that the 
multichannel images contain multiple regions of homogeneous texture. This is found to 
be the case in dealing with aerial photographs of natural terrain. 

The problem of multichannel image segmentation is addressed as an estimation 
problem for two regions of texture. Maximum likelihood {ML) and maximum a 
posteriori {MAP) region estimates using a Markov random field to model region 
transitions are developed. 

This chapter consists of two sections. The first section describes the linear 
filtering model and develops an expression for the image probability density function in 
terms of filter error residuals. The last section deals with the algorithm that is 
developed for the texture estimation and the multichannel image segmentation. The 
results of texture estimation and image segmentation are presented in Chapter IV. 

A. MODEL DEVELOPMENT 

In this Section, a 2-D multichannel autoregressive-moving average {ARMA) model 
[Ref. 5] is first discussed. Then, we concentrate on the multichannel autoregressive 
{AR) model with Gaussian white noise inputs [Ref 6]. The development parallels that 
in [Ref 1]. A multichannel image is represented by a vector signal (n,m) where 
(n,m) are spatial coordinates and the superscript h is an index representing the texture 
type. A 2-D multichannel image is shown in Figure 2.1 . A texture of type h is then 
modeled in general by a multichannel ARMA process defined by 
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( 2 . 1 ) 



P-1 Q-1 

(n,m) = “ X Z + Z it (n,m) 

i=0j=0 p 

(i,j) X (0,0) 



F\ (n,m) = t (n,m) + (2.2) 

for h = 0,1, n = 1,....,N, m = l,....,M, where and are set of Filter weighting 
coefficient matrices of size K by K. (n,m) are a set of independent identically 
distributed zero-mean random variables, £ is a constant representing the mean value of 
the image, and P is finite-extent mask covering the filtered points. F*’ (n,m), (n,m), 

and ^ are vectors of size K, the number of channels in the image. The matrices are 
the key parameters in the linear model. For a first quadrant filter with a P ^ Q region 
of support, there are PQ A^j^ matrices with A^qq = I, an identity matrix. 

For the auto-regressive {AR) or all-pole model we have and so that 

Equation 2. 1 reduces to 

P-1 Q-1 

t (n,m) = - I X A^.. t (n-i,m-i) + W (n,m) (2.3) 

i = 0j = 0 
(i,i) == (0.0) 

If the vectors f, w, and g represent an ordered set of the corresponding image points, 
then Equation 2.2 and Equation 2.3 are written in a matrix formulation as 

A(f - g) = W-Agfo (2.4) 

where A and Aq are matrices whose nonzero elements are derived from the terms Ajj in 
Equation 2.1 and fg represents a set of boundary conditions with support outside of the 
regions. Since the terms N (n,m) are independent with probability density function 
(PDF) p^^. One can solve Equation 2.4 for W and express the multivariate probability 
density function for the image conditioned on the boundary values as 
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Figure 2.1 2-D Multichannel Image Model. 
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(2.5) 



Pf| f,(f«fo)=y^|P,(A(f-g) + A^fo) 

= n f (n,m)) 

(n,m) e R 

where the notation E(n,m) is used to represent the ordered components of the vector 
A( f-g) + Aq Tq . If the boundary conditions, fg , are temporarily ignored, then 



E (n,m) = A (f - g) 



(2.6) 



and the terms E (n,m) in Equation 2.5 are computed from 



P-1 Q-1 

e'^ (n.m) = X E A^ij (n-i,m-j) (2.7) 

i=0 j = 0 
(i,i) (0,0) 

The filter of Equation 2.7 w'hich computes E ^ (n,m) from (n,m) is referred to as 
the 'prediction error fdter'. One can think of Equation 2.7 as producing an estimate or 
prediction F (n,m) of the image at point (n,m) and then forming an error E (n,m) as a 
difference F (n,m) — F (n,m) . This process is known as 2-D linear prediction and is 
fundamental to performing multichannel image segmentation. 

B. ALGORITHM DEVELOPMENT 

In the multichannel segmentation problem, it is assumed that the image consists 
of multiple connected regions of known texture types, but that the region boundaries 
and the number of regions are not known. The segmentation of the image is treated as 
a supervised learning problem, since the regions are considered to consist of known 
texture types. In this section, the multichannel image segmentation algorithm for 
textured images is discussed. 



12 



An overview of the method is as follows. Given a multichannel image of each 
texture, filter parameters are estimated by computing the covariance matrix from a set 
of data and solving the Normal equations corresponding to the model of Equation 2.3 . 
In this case the 'correlation method' of linear prediction is used to compute the 
covariance matrix. The filter parameters are derived from a statistical analysis of the 
textured images, because the image model discussed in the previous section is based 
upon statistical properties. 

Once the filter parameters are known the filters are used to perform the 
segmentation. The filter weighting coefficients are used to calculate the prediction 
errors E*’ (n,m) of two textures (n,m). Then, a maximum likelihood (ML) region 
estimate is developed using the prediction errors and the covariance matrices for image. 
The ML estimate is used as a basis to determine an approximate maximum a posteriori 
(MAP) region estimate. The MAP region estimation utilizes an underlying Markov 
structure for the region statistics to produce an accurate segmentation. 

1. Filter Parameter Estimation Method 

The prediction error filter is a finite-extent impulse response {FIR) or 
nonrecursive filter with selectable mask size and quart erplane region of support. It is 
always stable. The inverse of the filter used in Equation 2.3 is an all-pole filter (the AR 
filter). The filter parameter estimation problem requires calculating and A^ 

by statistical analysis of data in an estimation window containing the desired texture, 
where the quantity is a mean vector of the average gray level of the image in each 
of the channels, and is a covariance matrix of the multichannel white noise 



Z\=E[ (n,m) {W^ (n,m))T ] 



( 2 . 8 ) 



the term is also refered to as the prediction error covariance matrix, since in a 
linear prediction problem it represents the covariance of the quantities E (n,m) defined 
in Equation 2.7 . Since is not in general a diagonal matrix, we see that the 'white 
noise' is uncorrelated w’ithin each channel but correlated between channels. 

The covariance of the multichannel white noise (the prediction error 
covariance) and the filter coefficients will be obtained by estimating the correlation 
function of the image and solving a set of Normal equations as discussed below. The 
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correlation function itself is estimated from data in a window containing a sample of 
the desired texture. Two estimation windows are depicted in Figure 2.2 for two 
dilferent textures in the multichannel image. The reader should realize that the 
estimation windows do not have to come from the image to be segmented; they can be 
selected from any image containing the same type of texture. 



Figure 2.2 Typical estimation windows for two textures. 
a. Mean Vector Estimation 

In order to model the multichannel image by Equations 2.2 and 2.3 the 
mean vector of the multichannel image has to be estimated. Knowing M and selecting 
the stationary estimation windows of two desired textures of F^, the zero mean 2-D 
multichannel image, F, appearing in Equation 2.3 can be obtained by subtracting the 
mean vectors from the multichannel image. Thus, Equation 2.2 becomes 




F*’ (n,m) = F*’^ (n,m) - li' 



(2.9) 
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where M*‘ corresponds to g*' in liquation 2.2, and the term in liquation 2.9 arc 
computed from 




.M'-j , 



( 2 . 10 ) 



where 



1 X Y 

I ^^"2 '' 2 

^ 1 1 i\ (n,m) • (2. 1 1) 

n' m' n = Xjm=Yj 

is the mean estimate for the k**’ channel of the h* texture image. The limits Xj , 
X 2 , Yj , Y 2 represents the edges of the window which is of size N' by M'. Therefore 
n' = X 2 - Xj + 1 , m' = Y 2 - Yj + 1, andO Xj < X 2 ^ N, 0 ^ Yj < Y 2 
^ ,and h represents the two textures of the multichannel image. All variables used 
in Equation 2.11 are depicted in Figure 2.3 . 




Figure 2.3 Selection of the Size of the Windows for Mean Vector Estimation. 
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b. Correlation Function Estimation 

The correlation function of the zero mean, 2-D multichannel signal has to 
be calculated in order to estimate the multichannel white noise covariance or prediction 
error covariance, and the filter weighting coefficients, A*^jj . The theoretical 2-D 
matrix correlation function for lag " i,j " is given by 

R^'aj) = (R^* (-i, -j))T (2.12) 

= E [F^ (n,m) . ( F^ (n-i,m-j))^ ] 

and can be estimated from the multichannel signal by 

1 Vi2 m 2 

R^' (i,j) Y t (n,m) ( F*’ (n-i, m-j)T (2.13) 

N M n=nj m = mj 

where R^ (i,j) is a matrix of size K by K, and n^ , n 2 , m^ , m 2 are defined by 

0 ^ nj = max(Xj ,Xj + i ) < n 2 = min(X 2 , X 2 + i ) ^ N, and 
0 ^ mj = max(Yj , Yj + j) < m 2 = min(Y 2 , Y 2 + j ) ^ M. 

This matrix correlation function is used to form a larger block Toeplitz covariance 
matrix which is used to estimate the filter parameters. This is discussed next. 

c. Filter Coefficients and Prediction Error Covariance 

The prediction error filter weighting coefficients and the prediction error 
covariance must satisfy a set of linear equations known as the Normal Equations when 
the multichannel image is represented by the model in Equation 2.3 . 

Normal equations corresponding to Equation 2.3 can be written as 

[R] . [A] = [S] (2.15) 

where R is the correlation matrix for the signal, A is an appropriately ordered matrix 
of the filter coefficients, and S contains a single non-zero block which is the 
prediction error covariance. The matrix R has three levels of partitioning and for any 
rectangular region of support is block Teoplitz with block Toeplitz blocks. 

For a first quadrant filter with a P x Q region of support. The Normal 
equations that define Equation 2.15 have the specific form 
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•R(0) R(-l) . . 
R(l) R(0) . . 


. R(-P+l)' 
. R(-P + 2) 




< < ■ 

1 




1 

o O 
00 

1 


.R(P-l) R(P-2) . 


. . R(0) j 




AP * * , 




0 



(2.16) 



where 



R(i) = 



•R(i,0) R(i,-1) . . .R(i,-Q+1) 
R(i,l) R(i,0) . . .R(i,-Q + 2) 



LR(i,Q-l) R(i,Q-2) . . R(i,0) 

= (-i) 



(2.17) 



and where R(i,j) is the matrix correlation function described in Equations .2. 12 and 2.13 
The quantities A‘ and S® are defined by 



A' = 



A: 



i, 0 



A: 



i, 1 



A. 



i, Q - 1 



and 



(2.18) 



S = 




0 



( 2 . 20 ) 



17 



with Aqq= I and where Ajj and the partitions of S® are matrices of order K, the 
number of channels of image. 

2. Multichannel Segmentation Method 

An overview of the segmentation is as follows. By using a Gaussian 
probability density for the white noise in Equation 2.5 , one can develop explicit 
estimates for the density functions in terms of the prediction errors E (n.m) . From this 
one can form the conditions for ML and MAP estimation of the regions in the image. 
The theory leading to the estimates is explained below. A brief intuitive explanation of 
the process is given here. 

As mentioned earlier the prediction error filters in Equation 2.3 can be 
considered as predicting the intensity of a pixel in each channel from data in the region 
adjacent to the vector of pixels. The prediction errors ^ (n,m) are the outputs of the 
filters. In the segmentation algorithm the prediction error is normalized in an 
expression involving the corresponding prediction error covariance. These normalized 
errors are compared in an appropriate formula to obtain the ML region estimate. 
When an area of texture is processed by a filter that is not matched to the texture, the 
normalized prediction error can be expected to be high. ’When the same area is 
processed by the filter that is matched to the texture, the prediction error can be 
expected to be low. 

The 'maximum likelihood' region estimate, ML(n,m), of texture class is 
achieved based on the prediction error covariance and the error estimates of the two 
desired textures for each pixel in the multichannel image. The ML(n,m) region 
estimate assigns pixels to texture types without regard to the assignments of the 
adjacent pixels. Then, the 'maximum a posteriori' region estimate, MAP(n,m), of 
texture class is achieved for a pixel and a desired number of adjacent pixels of two 
textures of the ML region estimation result. The MAP region estimation uses the 
Markov model that refers to the above description. The form of the Markov model 
and ML and MAP region estimates are presented in detail below. 

a. Maximum Likelihood Region Estimation 

It is supposed that a multichannel image has many regions, but that each 
region contains only one or another of two texture types. Given these regions, one can 
w'rite the Equation 2.5 as 

p(f I Rj) = n p {E (n,m)) i = 1 , , q (2.21) 

hi 

(n,m) 
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where the is the probability density function for the white noise source of type hj 
hi 

within region Rj ,and q is the number of regions. When the white noise term is 
Gaussian with density function (mean 0 and covariance ) 

PU m exj, / - sT ‘ K ) (2.22) 

(2<t')' |j:„i 

then, taking minus twice the log of Equation 2.21 and applying Equation 2.22 , we 
obtain for an N by M pixel multichannel image 



- 2Jnp(£|R, ,Rj R^) 

= Z (lS‘h ("•m)!’'^ r ‘ IS‘h (n,m)l + ln|i:;,| + .. 

(n,m) e Rj 

+ Z ([^h (n.m)]T [Z\ ]- 1 [E\ (n.m)] + ln|i:\ |) 

(n,m) s Rj^ — NM ln2Tt 

(2.23) 



q . 

= I I (fs'h I-h )■ ‘ (n,m)l + ln|i;'^ I) 

i= 1 (n,m) € Rj 

- NM ln2it 



For maximum likelihood estimation, the number of regions q and the regions 
themselves are considered to be deterministic parameters of the density function. An 
ML estimate for these parameters is obtained by choosing values that maximize 
Equation 2.21 or, minimize Equation 2.23 . Since NMln27t is constant value, the 
Equation 2.23 is minimized if every point (n,m) in the multichannel image to a region 
Rj of type hj such that the term in brackets is minimum. Thus, one can WTite a ML 
region estimation for two textures as 



0 

MLj (n,m) > MLq (n,m) (2.24) 

1 



where 
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(n,m) = (n,m)]T [Z\ ]'^ [ 1 ^ (n,m)] + InlL^^I 



(2.25) 



for h = 0 , 1, where the number above or below the inequality indicates the region 
class to which the pixel (n,m) is assigned. When the class is 1, the ML(n,m) is assigned 
1 for a white pixel. Otherwise ML(n,m) is assigned 0 for a black pixel. Since the ML 
region estimate assigns pixels to black and white without regard to the assignments of 
adjacent pixels, this algorithm produces a number of false assignments and a somewhat 
'spotty' result. 

b. Maximum A Posteriori Region Estimation 

The 'maximum a posteriori' (MAP) region estimation utilizes the Markov 
model to describe the occurance of regions in the image. The combination of the linear 
filtering model with the Markov model results in an algorithm to achieve a MAP 
region estimation. For MAP estimation the regions are considered to be random 
quantities, and we maximize the probability for a given set of regions conditioned on 
our observation of the multichannel image. From Bayes rule, the a posteriori 
probability can be written as 



Pr [R, ,R 2 , I i ] 



P(£ |Ri .R 2 Rg ) • Pr[Ri . R. Rq ] 

p(f) 



(2.26) 



Since the denominator of Equation 2.26 does not depend -on the regions, maximum a 
posteriori (MAP) estimate for the regions can be obtained if the Rj are chosen to 
maximize the numerator 



Pil |Ri .R2 -Rq ) • Pr [Rj ,R2 ,Rq ] (2.27) 

One can define the " state " s(n,m) of point (n,m) as the region type to 
which that point has been assigned. In our development, the number of region types is 
assumed to be 2. Since the set of all possible state assignments for points in the image 
is one-to-one with the set of all possible divisions of the multichannel image into 
regions, the region estimation problem can be viewed as one of estimating the states of 
the points. It is assumed that the state of a point is stochastically dependent on some 
adjacent set of states m ^ symmetric support region, as shown in Figure 2.4 . 
Since S represents a chosen set of state assignments for all points in the multichannel 
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image, Pr [ S ] denotes the joint probability that the points in the image take on a 
chosen set of state assignments. The support set defines a neighborhood structure i.e. 
all elements in the support set are neighbors of each other. It can be shown that if the 
set of states S is a Markov random field, then the probability of S can be factored as a 
product of terms of functions depending only on the "cliques" of the support set . 
The cliques are defined as groups of points such that each set of points are neighbors 
of each other according to the support set. For a Markov random field, the probability 
of S can be written as 
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Figure 2.4 State support region of the point s for MAP estimation. 

Pr [ S ] = n Pr [ s(n,m) | S„ ^ ] (2.28) 

(n,m) 

where the terms in the product are additive functions defined on the cliques of 
and the product is over all cliques in S. One simple acceptable form for the terms in 
the product is 



Pr [ s(n,m) I ] = ^ exp [ s(n,m) {a + (t+t') + p^ (v + v') 

+ Yj (u + u') + Y 2 (2.29) 

where t = s(n-l,m), t' = s( n+ I, m), is the set of states as shown in Figure 2.4 , 
and D is a normalizing constant. One particular selection of the parameters, namely a 
= “ 4, Pj = P 2 = Yi “ Y 2 “ ^ leads to a particularly simple algorithm. In this 
case, we have 
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P*" [ * I ^n, m ] " exp ( X ( s( i, j ) - 1 / 2 )) 



(2.30) 




and 




(2.31) 



The second term in the numerator of Equation 2.26 can be replaced by 
Pr[S] . Thus maximizing Equation 2.26 is equivalent to minimizing 



One can define the MAP region estimate by combining Equations 2.22 , 
2.26 , 2.28 , and 2.32 as 



For Pr [ h I m J Equation 2.30 and Equation 2.31 computing the 

terms — 2 In Pr [hjS^^ ] is equivalent to counting the number of pixels in that 

have value 'h' and dividing by the total number of pixels in , and multiplying by 

an appropriate normalizing factor, KS.* A larger state support region is depicted 
in Figure 2.5 as an example. The side of the ^ must be an odd number . Although 
it does not necessarily lead to the true MAP estimate we find it convenient in practice 
to maximize Equation 2.33 term by term. That is, we require to use Equation 2.33 
without sum. Equation 2.33 can be solved by iteration using the maximum likelihood 



- 21n p( F I Rj ,....,Rq ) - 21n Pr [ S ] 



(2.32) 




0 

- 2 In Pr [ 1 I S„ „ ] < 




(2.33) 



+ ln|i:»„| - 21nPr[0|S„_„ 



Equat 

function. 
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state estimates obtained from liquation 2.24 as an initial set of states for MAP region 
estimate. 
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Figure 2.5 A Set of States I's, O's adjacent to S for MAP Region Estimate. 

The computational requirements of the MAP region estimation are reduced by storing 
the differences, 

MLD (n,m) - MLj (n,m) — MLQ(n,m) (2.34) 

incurred during the calculation of the ML region estimate. Substituting Equation 2.34 
in Equation 2.33 gives 

MLD(n,m) - 2 In Pr ( 1 | S„ ] + 2 In Pr [ 0 1 S„ ,„ ] < 0 (2.35) 

0 



At each iteration terms Pr [ h | S„ ] are evaluated based on the values of the states 
at the previous iteration. For our particular method of selecting Pr I h | S^^ , 

Equation 2.35 can be expressed as 



MLD(n,m) 

2 (number of state 



- KS 



1 pixels) — ( number of pixels in S^^ + I 



number of pixels in S 

^ n , ni 



1 

< 0 
> 

0 



(2.36) 
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There are two important points in Equation 2.36 in order to perform the maximum a 
posteriori image segmentation accurately. The value of the convergence factor, KS, 
must be assigned properly. If KS is assigned too small, the segmentation may not 
remove improperly classified pixels. On the other hand, if KS is assigned too large, 
correctly classified pixels could be changed. In addition, the size of S^^ must be large 
enough. Otherwise, false assignments produced by the initial ML segmentation may 
not be removed. 
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III. KARHUNEN-LOEVE TRANSFORMATION AND ONE-CHANNEL 

IMAGE SEGMENTATION 

In this chapter, the models and relevant algorithms are presented to perform the 
Karhunen-Loeve (K-L) transformation [Ref 3] and one-channel image segmentation 
utilizing the techniques of linear prediction [Ref I]. Since linear prediction techniques 
are presented in detail for multichannel image segmentation in the previous chapter, 
this discussion concentrates on the K-L transformation. The K-L model and algorithm 
are first developed to reduce the three-channel color problem to a one-channel 
problem. Then, a one-channel segmentation procedure is presented that is based on the 
same model previously discussed. The results of the K-L transformed one-channel 
image segmentation are presented and compared with the multichannel image 
segmentation in Chapter IV. 

A. MODEL DEVELOPMENT 

The K-L transformation developed in this section is based on the statistical 
properties of an image. This transformation provides an energy compaction between 
channels of a color image. That is, most of the color image energy is compacted into 
one channel, and the transformed image channels are uncorrelated. If the multichannel 
image and transformed multichannel irhage are expressed in vector form. The K-L 
transformation is given by [Ref 7] 

2(n,m) =[A] F{n,m) (3.1) 

where F ( n , m ) is the original multichannel image, g ( n , m ) is the transformed 
multichannel image, and A is the K-L transformation matrix, whose rows are 
eigenvectors of the between-channel correlation matrix R defined by 

R = E [ F ( n , m ) ( n , m ) ]. (3.2) 

The between-channel correlation matrix of the transformed image is 

A= E[g(n,m)2'^(n,m)] (3.3) 

= [ A][R][ A]-l 
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where the matrix A is a diagonal matrix of the form 



A = 



■ Xj 0 0 

0 0 
. 0 0 X3 



(3.4) 



and the represent eigenvalues of the between-channel correlation matrix. The 
eigenvalues are ordered such that 



Xj > X 2 ^ X 3 ^ 0 (3.5) 

The importance of this property is that each eigenvalue Xj^ is equal to the variance of 
the k * channel of the transformed multichannel image whose channels • are 
uncorrelated. Then, it is a well knoum property of multivariate statistics that the total 
variability of the color image has the form [Ref 3] 

3 

Xj = X X^ (3.6) 

k= 1 

which relates the total variability to the decorrelated component variations, Xj^ . One 
can observe that often the Xj^ values have a wide range of magnitudes, and the first 
component, Xj , will be sufficient to approximate X-p with only a small percentage of 
error. This becomes the key idea for the use of the K-L transformation. Table 1 shows 
the energy distribution between the transformed color image channels of two test 
images. 

Indeed, in a 3-channel image one can often find that the first channel of the K-L 
transformed image is sufficient to account for 99 percent of all the variability. That is, 
it is typical to find that 

Xj = 0.99 Xj (3.7) 

The Karhunen-LoSve transformation provides the best energy compaction [Ref 8] and 
the advantage of this transformation is computational savings. We will see later that 
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TABLE 1 



ENERGY DISTRIBUTION BETWEEN TRANSFORMED IMAGE 

CHANNELS 





Percentage of Energy 


in Channels 




Qi Q2 


O' 


Image I 


99.13 0.84 


0.03 


Image 2 


99.15 0.78 


0.07 



one-channel image segmentation achieved by processing only the first- channel of the 
K-L transformed color image will be very close to the result obtained on the original 
image with a multichannel algorithm, but will be obtained at only one ninth of the 
computational cost. 

B. ALGORITHM DEVELOPMENT 

In this section, the algorithm to perform K-L transformation of color images is 
presented. Since the K-L transformation is based upon the between-channel correlation 
matrix, R , of the color image , the between-channel correlation matrix is first 
determined. Then the eigenvectors of the between-channel correlation matrix are 
computed. Finally, the K-L transformation is formulated using the transpose of the 
eigenvector matrix. 

1. Correlation Function Estimation 

In order to determine the K-L transformation of Equation 3.1 , the between- 
channel correlation matrix of the color image must first be estimated. Our estimate for 
the between-channel correlation matrix is given by 



1 N-1 N-1 

R = m F ( n , m ) ( F ( n , m 

N^ n = 0 m = 0 



(3.8) 
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where the image size is N by N, and the between-channel correlation matrix size is K 
by K. 

2. Karhunen-Loeve Transformation Matrix 

Since the rows of the K-L transformation matrix are the eigenvectors, E, of 
the between-channel correlation matrix, the eigenvectors must be calculated from the 
between-channel correlation matrix of the color image. The K-L transformation 
matrix is then obtained using the transpose of the eigenvector matrix. That is 

(3.9) 




3. Karhunen-Loeve Transformation 

Since the K-L transformation matrix satisfies Equation 3.3 , then the K-L 

transformation is represented by Equation 3.1 with A is given by Equation 3.9, where 

T 

e j , i = 1,2,3 are the eigenvectors of R. More explicitly, the components of g are 
determined from the components of F at any pixel ( n , m ) as 



■Qj ( n , m ) 






'F j ( n , m ) 


Q2 ( n , m ) 


= 


^ ^ ^ 


F2 ( n , m ) 


Q3 ( n , m ) 




^ ^ 


F3 ( n , m ) 



(3.10) 



C. ONE-CHANNEL IMAGE SEGMENTATION 

The one-channel segmentation procedures presented in this section are based on 
the same model that w'as described in detail for the multichannel image in Chapter II. 
A single channel image is used instead of a three channel image. That is, K is always 
equal to 1 in the Equations of Chapter II. As a result the vector and matrix quantities 
become scalars and the equations are simplified. 

In order to model the single image by Equations 2.2 and 2.3, the mean of the 
image is estimated using the Equation 2.11 . Then the correlation function is estimated 
in the same manner as in Section II-B.b where R(i.j) is a scalar instead of a matrix. 
This procedure is followed by estimation of the filter coefficients and the prediction 
error covariance. The equations in Section II-B.c are then used to determine the filter 
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coeHicients, Ajj , and the the prediction error covariance, , now a scalar value. 
Finally, the one-channel segmentation algorithm is applied. The method is the same as 
that described in detail for the color images in Section II-B.2. However, instead of 
Equation 2.24 and 2.33 the foil owning simplified relations are used for the maximum 
likelihood and the maximum a posteriori region estimates. A maximum likelihood 
region estimate for a single channel image of the pixel is given by [Ref 1] 



(n,m))^ 



t;1 



W 




(3.11) 



and the maximum a posteriori region estimate is given by 




■“ w 



n, m 



(3.12) 



where E^ and E® are the result of applying the linear predictive filters to the first 
channel of the K-L test image. 
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IV. RESULTS AND COMPARISON OF THE METHODS 



In this chapter, the results of the multichannel image segmentation, the 
Karhunen-Loeve transformed one-channel image segmentation, and the K-L 
transformed 3-channel image segmentation are presented and compared. 

The digitized image size used in this work is 128 by 128 pixels with gray levels 
represented on a scale of 0 to 255 (8 bits). A digitized color photograph of a rural area 
containing trees (the green region) and fields (the yellow region) are shown in Figure 
4.1 . A quarter-plane filter for each texture class (2x2 pixels) w'as designed and 
applied to the color image to achieve the multichannel image segmentation. A state 
support region of 7 by 7 pixels w'as used for MAP region estimation. The results of the 
maximum likelihood (ML) and the maximum a posteriori (MAP) segmentations of 
Figure 4.1 image are shown in Figures 4.2 and 4.3 respectively. The segmentation 
results show the field regions as black and tree regions coded as w'hite. The ML result 
(Fig 4.2) is spotty, but the true tree and field regions are distinguishable. The MAP 
segmentation result (Fig 4.3) of Figure 4.1 image is quite clear. The MAP 
segmentation w'as not able to remove a few improperly classified points in the left side 
of the field region. Since there is an inherent ambiguity in the estimation of the region 
boundaries due to the finite size of the masks, the edges are expectedly somew’hat 
rough. Figure 4.4 shows another color image of a rural area containing treds and 
fields. Figures 4.5 and 4.6 show the results of the ML and the MAP segmentation of 
the Figure 4.4 image. The ML estimation was achieved using the filters designed for 
the image of Figure 4.1 . The result of the ML segmentation (Fig 4.5) is again quite 
spotty, i.e. Although two regions are discernible. The MAP algorithms segmented the 
regions quite accurately as shown in Figure 4.6 . The MAP estimates presented above 
converged after 10 iterations and KS was assigned to 100. 

In the results presented above the ML procedure produced a poor result with a 
lot of false regions. This is due to the lack of prior information about region 
connectivity. On the other hand, MAP estimation using the Markov model to represent 
region transitions produced results that w’as quite accurate. 

Figure 4.7 shows the first channel of the Karhunen-Loeve (K-L) transformed 
image of the Figure 4.1 . The image size is 128 by 128 pixels with the scaled gray levels 
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within the intensity range of the display. The result of the one-channel ML 
segmentation of Figure 4.7 is depicted in Figure 4.8 . The one-channel ML 
segmentation result is quite spotty, but the two regions are perceptually discernible. On 
the other hand, the MAP segmentation result (Fig 4.9) of Figure 4.7 is quite smooth, 
although there are a few incorrectly classified points in both regions. 

The first channel of the Karhunen-Loeve transformed image of Figure 4.4 is 
shown in Figure 4.10 , and Figures 4.11 and 4.12 show' the results of the one-channel 
ML and MAP segmentations of Figure 4.10 respectively. The one-channel ML 
segmentation result (Fig 4.11) has a lot of spots, but both segmented regions are 
distinguishable. The maximum a posteriori segmentation (Fig 4.12) of Figure 4.10 is 
quite smooth, but again there are a few incorrectly classified sets of pixels in both 
regions. 

The results of multichannel image segmentation and the Karhunen-Loeve 
transformed one-channel image segmentations were presented consecutively. These 
results show that the ML estimation of the Karhunen-Loeve transformed single 
channel image is much more spotty than the ML estimation of the multichannel image. 
But, the MAP estimation results of the K-L transformed single channel images are as 
clear as the MAP segmentation results of the multichannel images. 

In summary the results presented in this chapter show that the best segmentation 
result is provided by the multichannel image segmentation method. However the 
segmentation results of the Karhunen-Loeve transformed single channel images are 
very close to multichannel image segmentation results, i.e. Because most of the color 
image energy (99 percent) is compacted into the first channel. 
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Figure 4.1 Two-Texture Color Image. 
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Figure 4.2 ML Region Estimation of Figure 4.1 Image. 




Figure 4.3 MAP Region Estimation of Figure 4.1 Image. 




( 

Figure 4.4 Color Image Containing Two-Texture. 
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Figure 4.5 ML Region Estimation of Figure 4.4 Image. 




Figure 4.6 MAP Region Estimation of Figure 4.4 Image. 
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Figure 4.7 First Channel of K-L Transformation of Figure 4.1 Image. 
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Figure 4.8 ML Region Estimation of Figure 4.7. 
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Figure 4.10 K-L Transformed Single Channel Image of Figure 4.4. 
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Figure 4.1 1 ML Segmentation of Figure 4.10. 




Figure 4.12 MAP Segmentation of Figure 4.10. 
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V. CONCLUSIONS 



The segmentation of terrain images is an important part of image analysis 
methods for military and civilian applications. The work in this thesis utilized a 2-D 
stochastic linear filtering model and compared algorithms for multichannel image 
segmentation of color images. Two levels of structure were used for multichannel 
segmentation development. The fundamental structure based on the linear filtering 
concepts represents the texture in local regions of terrain. Superimposed on this 
structure is a Markov random field that describes transitions from one region type to 
another. The segmentation W'as considered as a region estimation problem and 
maximum likelihood and maximum a posteriori region estimatation methods were 
developed. The ML region estimation produced a spotty result, but the MAP region 
estimation produced quite accurate results for the multichannel and single channel 
image. 

The other piece of work developed in this thesis was the Karhunen-Loeve 
transformation model that based on the statistical characteristics of color image.The 
one-channel image segmentation was then applied to the first channel of the 
Karhunen-Loeve transformed color image to see the effectiveness of the K-L 
transformation for segmentation. 

We observed that multichannel image segmentation results were quite accurate. 
Similarly the results of the K-L transformed one-channel image segmentation were very 
smooth. In summarv’ the results of both segmentation methods were veiy^ close to each 
other, and the K-L transformation is very effective for segmentation. 
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APPENDIX A 
RELAXATION METHOD 

The Relaxation Method algorithm utilizes a set of iterative numerical techniques 
to compute a posterior probability of the pixel (n,m) from a prior probability of the 
same pixel and a set of prior probabilities of the adjacent pixels. The prior probabilities 
are estimated from a 2-D image data using the linear filtering model (see Equation 
2.25). 

The Relaxation formula is defined [Ref 9] by 
k + 1 K" Pn‘‘ (n,m) 

Pjj ^ Mn,m) = Ap( — ) (A.l) 



where k is the number of the iteration, Xy® is the relaxation factor, Pj|*^ (n,m) are a set 
of prior probabilities, T is the number of textures, and S is the number of pixels. The 
updated estimates Pj.^ ^ (n,m) are obtained by averaging all of the terms in 
parentheses. The relaxation factor, X.j® , in Equation A.l is given by 

V = J c(i,j|s,t)P K{n,m) (A.2) 

*1 t = 1 



where c(i,j|s,t) is a nonnegative compatibility function whose value is small if the 
neighboring pixel is black when the estimated pixel is white, otherwise its value will be 
large. Xjj* is used to update the probability Pjj^ (n,m). Note that Xjj® is large if the 
compatibilities c(i,jls,t) are large and the probabilities P^^^ are high, otherwise Xjj* will 
be small. 

The results of the Relaxation Method segmentation are shown in Figures A.l 
and A.2 . The Figure A.l presents the segmentation of Figure 4.7 with the 
compatibility c(i,x|s,x) is equal to 0.1, where x can be 0 or 1. The result is ver>’ spotty 
but both regions are discernable. The Figure A.2 gives the segmentation with the 
compatibility c(i,x|s,x) is equal to 0.9. This result is better than the previous result, but 
there are still several spots especially in the field region. Both results are obtained after 
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lU iterations. I'igure 4.9 shows the result of MAP segmentation of f igure 4.7 . The 
MAP segmentation result is much more accurate than the Relaxation method 
segmentation. 
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Figure A. 2 Relaxation Method Result with c = 
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APPENDIX B 

FILTER PARAMETER^ ESTIMATION AND MULTICHANNEL 

SEGMENTATION 



The interactive program, FLTRl, estimates two sets of filter parameters from a 
2-D three-channel image. The corresponding algorithms were presented in Section 
II.B.l. In order to run this program, the user must input the required program 
parameters in the order listed below 

* The number of filter rows, P. Maximum is 4. 

* The number of filter columns, Q. Maximum is 4: 

* The number of rows, N, in each image channel. Maximum is 12S. 

* The number of columns, M, in each image channel. Maximum is 128. 

* The number of channels, K, in image. 

* The filenames of the image channels. 

* The coordinates of the esumation windows of textures. 

* Tw'o output filenames for the estimated filter parameters. 

The mean vectors of the data in the two estimation window’s is first computed by 
FLTRl. Then, the program subtracts the mean from the image and determines the 
correlation functions. After calculating the correlation matrix, the program determines 
the covariance matrix, , using the equation 



[R][B] = [Si] (B.l) 

where B is a dummy matrix that has the same dimension with A matrix, Bgg = [ ]‘ ^ 

and Sj is all zero except for an identity matrix in its first partition. The covariance 
matrix must satisfy the relation 



B, 



00 



I. = I 

w 



(B.2) 



Finally, the filter weighting coefficients are estimated using 



[A]= [B].[i:] 



(B.3) 



44 



MULTICHANNEL IMAGE SEGMENTATION PROGRAM 



The interactive program, SGMTl, segments a color image with two textures 
when given the three-channel image, the image dimension, and two sets of filter 
parameters. Again, the user must input the required program parameters to run the 
program. These parameters have to be in the order listed below : 

* The number of rows, N, in the image channels. 

* The number of columns, M, in the Tmage channels. 

* The number of filter rows, P. 

* The number of filter columns, Q. 

* The number of textures, T. in the image. 

* The number of channels, K, in the image. 

* The filenames of the image channels. 

* The filenames of the filter parameters. 

* The output filename of the ML segmentation. 

* The output filename of the MAP segmentation. 

This program estimates the errors of two textures using Equation 2.7, then performs 
the ML segmentation and the MAP segmentation using Equations- 2.25 and 2.33 . The 
convergence factor KS, and the size of must be assigned properly by user to 

perform the MAP segmentation accurately. 
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APPENDIX C 

KARHUNEN-LOEVE TRANSFORMATION 



The interactive program, KRLV, implements the Karhunen-Loeve transformation 
algorithms described in Section III.B. This program requires the user to assign the 
program parameters in the order listed below : 

* The number of channels, K, in the image. 

* The number of rows and columns, N, in each channel. 

* The filenames of the image channels. 

* The output filenames of the transformed color channels. 

The correlation matrix is first calculated. Then the program calls the IMSL 
subroutine EIGRS to calculate the eigenvalues and the eigenvectors of the correlation 
matrix. Finally, the transformed color image is obtained by multiplying the transpose 
of the eigenvectors matrix by the original color image. The program scales the 
transformed color image for display on the COMTAL image prossessing system. 



ONE-CHANNEL IMAGE SEGMENTATION PROGRAM 

The program, FLTR2, calculates two sets of filter parameters from a 2-D single 
channel image. The user must input the required program parameters to run the 
program. These parameters except the number of channel, K, are given in Appendix B. 

The program, SGMT2, implements a single channel image segmentation. The 
required program parameters given in Appendix B have to be assigned to run the 
program. Equations 3.11 and 3.12 are used to perform the ML region estimation and 
the MAP region estimation. 
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APPENDIX D 



COMPUTER PROGRAMS FOR MULTICHANNEL IMAGE 

SEGMENTATION 



********************************************************************** 



* 
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PROGRAM FLTRl 

PURPOSE To develop two sets of filter parameters from a 2-D 
three-channel input image. These parameters are the 
mean vector, the covariance, and the filter coeffici- 
ents of the image. 

REQUIRED TMSL ROUTINES 

LEQT2F, LUDATN, LUELMN, LUREFN 

IMPLEMENTED BY LTJG TIMUR KUPELI Nov 1986 
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**** VARIABLE DEFINITIONS **** 






Input 
Input 
Rows , Co 
Rows , Co 






Image data in byte format, 
image data in real*4 format, 
umns of the linear filter.- 
umns of the image data. 

The number *of channels in image. 

The number of filters for processing. 

’’ Input ] Filename of the image channels. 

Output ] Filename of the filter parameters. 
Array of mean vectors of estimation windows. 
The correlation matrix. 

Identity matrix. • 

The covariance matrix 

The filter coefficients matrix. 

The estimation window size. 

Row, column delay(shift) of correlation. 



B INPUT 
F INPUT 

N ' 8 
K 
T 

IMAGE 
FILTER 
MEAN 
R 

IMATRX 
KW 
A 

ISIZE 
NO, MO 

This program uses 2 by 2 pixels filter. If user wants to use 
different size filter , the dimensions of the following 
variables must be modified. 

R- SI, A, WAREA 

For exemple , for 4 by 4 pixels filter, the dimensions must be 
R(48,48), SI(48,3), A(48,3), WKAREA(2448) 



INTEGER*4 P , Q , ROW , COL , STARTN ( 2 ) , STARTM ( 2 ) , ENDN ( 2 ) , ENDM ( 2 ) , 

* RNROW,RNCOL,RROW,CCOL,X,Y,K,IER,QK,PKQ,L,J,JJ,JJJ, 

* HNNO , HMMO , HCOL , HROW , LNNO , LNMO , LCOL , LROW , RN , RM , T , 

* NO , MO , RRN , RRM , ROWl , COLl , MMO ,NN0 , IDGT , I 

REAL*4 SUM,TEMP,FINPUT(128,128,3) ,MEAN(2,3) ,WKAREA(180) , 

* R(12,12) ,SI(12,3) ,IMATRX(3,3) ,B00(3,3) ,SUM1 ,SUM2,SUM3, 

* KW(3,3) ,A(12,3) ,ISIZE 

CHARACTER*50 IMAGE (3 ), FILTER (2) 

BYTE BINPUT(128) 



GET PROGRAM INPUT PARAMETERS. 
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noon o o ooo o ooo 



TYPE 10 

10 F0RMAT(' ENTER NUMBER OF FILTER ROWS FROM 2-4 DESIRED s',$) 

READ 11, P 

11 F0RMAT(I3) 

TYPE 12 

12 FORMATC enter number of filter columns from 2-4 DESIRED 

READ 11, Q 

TYPE 13 

13 FORMATC ENTER NUMBER OF ROWS IN IMAGE ' ,$) 

READ 11, N 

TYPE 14 

14 FORMATC ENTER NUMBER OF COLUMNS IN IMAGE ',$) 

READ 11, M 

TYPE 15 

15 FORMAT?' enter NUMBER OF CHANNELS [ MAX = 3 ] IN IMAGE ;',$) 
READ(*,11) K 

GET THE MULTICHANNEL IMAGE 
DO 16 J = 1 , K 
WRITE (*,17) J 

17 FORMATC ENTER FILENAME OF IMAGE CHANNEL ' 13, $) 

READ(*,18) IMAGE(J) 

18 FORMAT (A50) 

CONVERT THE IMAGE FROM BYTE FORMAT- TO REAL*4 FORMAT 

OPEN(UNIT=l,FILE=IMAGE(J),STATUS='OLD' , ACCESS = 'DIRECT') 

DO 180 ROW = 1 , N 

READ (I'ROW) (BINPUT(COL) ,COL=l,M) 

DO 181 COL = 1 , M 

TEMP = BINPUT(COL) 

IF (TEMP. LT. 0.0) THEN 
TEMP = TEMP +256 
END IF 

F INPUT (ROW, COL, J) = TEMP 

181 CONTINUE 

180 CONTINUE 

CLOSE (UNIT=1) 

16 CONTINUE 

GET 'T' AREAS FOR WHICH FILTERS ARE DESIRED AND OUTPUT FILENAMES 
FOR EACH AREA'S FILTER COEFFICIENTS AND COVARIANCE MATRIX. 

DO 19 I = 1,T 
WRITE (*,20) I 

20 FORMATC ENTER UPPER-LEFT ROW FOR AREA ',12,':',$) 

READ(*,11) STARTN(I) 

WRITE (*,21) I 

21 FORMATC ENTER UPPER-LEFT COLUMN FOR AREA' ,12, ' ; ' ,$) 

READ(*,11) STARTM(I) 

WRITE(*,22) I 

22 FORMATC ENTER LOWER-RIGHT ROW FOR AREA ' , 12 , ' : ' , $) 

READ(*,11) ENDN(I) 

WRITE (*,23) I 

23 FORMATC ENTER LOWER-RIGTH COLUMN FOR AREA ' , 12 , ' : ' , $) 

READ(*,11) ENDM(I) 
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24 



C 

C 

C 

C 



26 

25 



C 

C 

C 



Z1 



29 

28 

271 

C 

C 

C 

C 

C 



291 



WRITE (*,24) I 

FORMATC enter output file-spec for FILTER' ,12, 

READ(*,18) FILTER(I) 

FIND THE MEAN VECTOR OF ESTIMATION WINDOW , T, AREAS OF IMAGE 
CHANNELS. THE MEAN VECTOR CONSISTS OF THE MEAN FOUND EACH CHANNEL 



ISIZE = (ENDN(I) - STARTN(I) + 1)*(ENDM(I) 
SUMl =0.0 
SUM2 =0.0 



SUMS =0.0 

DO 25 L = STARTN(I) ,ENDN(I) 

DO 26 J = STARTM(I),ENDM(I) 
SUMl = SUMl + FINPUT(L,J,1 
SUM2 = SUM2 ' ' 

SUMS = SUMS 
CONTINUE 
CONTINUE 



+ FINPUT(L,J,2 
+ FINPUT(L,J,S; 



MEAN(I,1 

MEAN(I,2 

MEAN(I,S 



= SUMl / ISIZE 
= SUM2 / ISIZE 
= SUMS / ISIZE 



WRITE(^,27) (MEAN(I,II),II=1,K) 
FORMAT(' MEAN: ' ,SF9.S,$) 



STARTM(I) + 1) 



• CORRECT THE IMAGE TO BE ZERO MEAN 



DO 271 J = 1 , K 

DO 28 L = STARTN(I), ENDN(I) 

DO 29 LL = STARTM(I) ,ENDM(I) 

FINPUT(L,LL,J) = FINPUT(L,LL,J) - MEAN(I,J) 

CONTINUE 

CONTINUE 

CONTINUE 

DETERMINE THE 2-D CORRELATION FUNCTION OF THE IMAGE CHANNEL. 

THE CORRELATION MATRIX APPROPRIATE TO THE INPUT INPUT PARAMETERS 
OF P, K, Q 

WRITE(*,291) I 

FORMATC CORRELATION. MATRIX' ,I2,$) 

PKQ = P * K * Q 
QK = Q * K . 

DO SO RNROW = 1,P 

X = QK * (RNROW-1) 

DO SI RNCOL = 1,P 

NO = RNROW - RNCOL 
y = QK * (RNCOL-1) 

DO S2 RMROW = 1,Q 

RN = K * (RMROW - 1) + X 
DO SS RMCOL = 1,Q 

MO = RMROW - RMCOL 

RM = K * (RMCOL - 1 ) + Y 



C LNNO = STARTN(I) + NO 

LMMO = STARTM(I) + MO 
HNNO = ENDN(I) + NO 
HMMO = ENDM(I) + MO 
C 

LCOL =MAX0 ( STARTM ( I ) , LMMO ) 
HCOL =MIN0(ENDM(I),HMM0) 
LROW =MAXO ( STARTN ( I ) , LNNO ) 
HROW =MIN0(ENDN(I),HNN0) 

C 

DO ISS R0W1= 1 , S 
RRN = RN + ROWl 

DO 23S COLl = 1 , S 
RRM = RM + COLl 
SUM = 0.0 

DO SSS ROW = LROW , HROW 
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433 

333 

233 

133 

33 

32 

31 

30 



NNO = ROW - NO 

IF (NNO .LT. LROW) THEN 

RROW = MAXO(ROW,NNO) 
ELSE IF (NNO ,GT. HROW) THEN 
RROW = MIN (ROW, NNO) 



331 

332 
330 



ELSE 
END IF 
DO 433 



RROW = NNO 



COL = LCOL , HCOL 
MMO = COL - MO 
IF (MMO .LT. LCOL) THEN 
CCOL = MAXO(MMO,COL) 

ELSE IF (MMO .GT. HCOL) THEN 
CCOL = MINO(MMO,COL) 

ELSE 

CCOL = MMO 

END IF 

SUM = SUM + F INPUT (ROW, COL, ROWl) * F INPUTS AVEW, CCOL, COLl) 

CONTINUE 

CONTINUE 

R(RRN,RRM) = SUM / ISIZE 
CONTINUE 
CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

THE FOLLOWING FORMAT MUST BE MODIFIED TO USE DIFFERENT FILTER SIZE 
THEN 2 by 2 PIXELS. 

DO 330 L = 1 , PKQ 
WRITE?*, 331) (R(L,J),J=1,PKQ) 

WRITE?*, 332) 

FORMAT?' ’,12F6.0) 

FORMAT? ' ' ) 

CONTINUE 



RESET THE IDENTITY MATRIX. 



DO 36 J = 1 , K 

DO 37 L = 1 , K 

IF (J.EQ.L) THEN 

IMATRX(J,L) =1.0 
ELSE 

IMATRX(J,L) =0.0 
END IF 

37 CONTINUE 

36 CONTINUE 

RESET SI(J,L) TO HAVE [ I ] IN FIRST PARTITION AND [ 0 ] IN ALL OTHERS 

DO 38 J = 1 , PKQ 
DO 39 L = 1 , K 

IF (J.EQ.L) THEN 
SI(J,L) = 1.0 
ELSE 

SI(J,L) = 0.0 
END IF 

39 CONTINUE 

38 CONTINUE 

SOLVE EQUATION [ R ] * [ B ] = [ SI 1 . NOTE THAT THE BELOW 
CALL TO IMSL ROUTINE LEQT2F, [ B ] RETURNS IN [ SI ] 
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IDG = 3 
C 

CALL LEQT2F (R ,K,PKQ ,PKQ , SI , IDG,WKAREA, lER) 

SOLVE EQUATION [ BOO ] * [ KW ] = [ I ] FOR COVARIANCE 
MATRIX [ KW 1 AFTER COLLING IMSL ROUTINE LEQT2F, [ KW ] 
RETURNS IN [ IMATRX ]. 

DO 45 J = 1 , K 

DO 46 JJ = 1 , K 

B00(J,JJ) = SI(J,JJ) 

46 CONTINUE 
45 CONTINUE 

CALL LEQT2F ( BOO , K , K , K , IMATRX , IDG , WKAREA , lER) 

SOLVE FOR THE FILTER COEFFICIENTS , [ A ] = [ B ] * [ KW ] , 
WHICH IN THE PROGRAM IS [ A ] = [ SI ] * [ IMATRX ] 

DO 47 J = 1 , PKQ 
DO 48 JJ = 1 , K 
TEMP =0.0 
DO 49 JJJ = 1 , K 

TEMP = TEMP + SI(J,JJJ) * IMATRX(JJJ,JJ) 

49 CONTINUE 

A(J,JJ) = TEMP 

48 CONTINUE 

47 CONTINUE 
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WRITE ( * , 49 1 ) ( ( IMATRX ( J , J J ) , J J=1 , K ) , J=1 , K ) 
FORMATC ■ ,<K>(F6.2,3X)) 

WRITE(*,53) I 

FORMATC FILTER COEFFICIENTS 12, $) 
WRITE(*,491) ( (A(J,JJ) ,JJ=1,K) ,J=1,PKQ) 



WRITE OUT MEAN, COVARIANCE MATRIX, AND FILTER COEFFICIENTS TO THE 
USER INPUT FILE. 



* 



C 



495 



C 

C 

19 
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OPEN (UNIT=2,FILE=FILTER(I) ,STATUS='NEW' , CARRIAGECONTROL= ' LIST ' , 
FORM= ' FORMATTED ' ) 

WRITE(2,495) (MEAN(I , J) , J=1 ,K) 

FORMAT (FIO. 4) 

WRITE (2, 495) ( ( IMATRX( J , JJ) , JJ=1 ,K) , J=1 ,K) 

WRITE(2,495) ( (A( J , JJ) , JJ=1 ,K) , J=1 ,PKQ) 

CLOSE (UNIT=2) 

CONTINUE 

WRITE(*,55) 

FORMATC THE PROGRAM FLTRl IS OVER',$) 

STOP 

END 
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PROGRAM SGMTl 

PURPOSE To segment 
the image, 
sions, and 
be used. 



a color image with two textures given 
the image dimensions, the filter dimen- 
the two sets of filter parameters to 



REQUIRED IMSL ROUTINES 

LINV3F, LUDATN, LUELMN 

IMPLEMENTED BY LTJG TIMUR KUPELI Nov 1986 



* 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 
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**** VARIABLE DEFINITIONS **** 



BINPUT 

ML 

MAP 

FNAME 

IMAGE 



K 

KW 

MEAN 

ERROR 

A 

TEXTURE 

COl 

CIO 

IN , IM 
PI , Q1 
TMAX 
IK 



Input ] Image data in the byte format. 

Output 1 The result of ML segmentation in byte format. 
Output j The result of MAP segmentation in byte format. 
Input ] Filename of filter parameters set. 

Filename of the image channel. 

Rows , Columns of filter. 

Rows , Columns of image. 

Number of channel of image. 

[ Input ] The covariance matrix. 

[ Input J Mean vector of estimation windows. 

Prediction error estimation 
[ Input ] Filter coefficients matrix. 

Zero-mean image data in real*4 format. 

Number of removed false points in the first texture. 
Number of removed false points in the other texture. 
Maximum number of rows and columns in the image. 

Maximum number of rows and columns in the filter. 

Maximum number of filters for processing. 

Maximum number of channels in image. 



INTEGER IN,N,IM,P1,P,Q1,Q,TMAX,T,PQ,R0W,C0L,I,J,JJ,JJJ,L,LL,KK, 

* LLL,LLLL,COUNT,LI,HI,LJ,HJ,K,PP,QQ,IK,II,III,C01,C10,M 

REAL KW(1 :3,1;3,1;2),MEAN(1:3,1:2) ,TEMP,SUM1,SUM2,PML11,LN(2) , 

* ERR0R(1 : 128,1 : 128,1 : 3,1:2 ) ,PML1,PML2,PML(1 : 128, 1:128), 

* AREA,KS,A(1:3,1 :3,1:2,1:2,1:2) ,D1,D2,KW1(1:3,1:3) ,KW2(3,3) , 

* EKWl (1,1:3) ,EKW2(1,1:3),PML22,DD2,TEXTUR(1:128,1:128,3,3) 

* WKAREA(6KAA(1:128, 1:128) 

CHARACTER*50 IMAGE (1:3), FNAME (1:2) ,MLTEST , MAPTEST 

BYTE BINPUTd : 128, 1:128, 1 :3) ,ML(1 :128, 1 :128) ,MAP(1 : 128 , 1 : 128) 

* MLI(1:128, 1:128) 

IN = 128 
IM = 128 
PI = 4 
Q1 = 4 
TMAX = 2 
IK = 3 

GET THE INPUT PARAMETERS OF THE PROGRAM 

1 WRITE(*,2) IN 

2 FORMATC ENTER THE NUMBER OF ROWS IN IMAGE. LIMIT OF ', 13 ,':',$ ) 

READ(*,3) N 

3 FORMAT (13) 

IF((N.LT.l) .OR. (N.GT.IN)) GOTO 1 



non on on non 



4 WRITE(*,5) IM 

5 FORMATC ENTER THE NUMBER OF COLUMNS IN IMAGE. LIMIT OF ' , 13 , ' : ' , 

READ(*,3) M 

IF((M.LT.l) .OR. (M.GT.IM)) GOTO 4 
C 

6 WRITE(*,7) PI 

7 FORMATC ENTER THE NUMBER OF ROWS IN FILTER. LIMIT OF',I3, ■;',$) 

READ(*,3) P 

IF((P.LT.2) .OR. (P.GT.Pl)) GOTO 6 
C 

8 WRITE(*,9) Q1 

9 FORMATC ENTER THE NUMBER OF COLUMNS IN FILTER. LIMIT OF ' , 13 , ' : ' , 

READ(*,3) 0 

IF((Q.LT.2) .OR. (Q.GT.Ql)) GOTO 8 
C 

10 WRITE (*,11) TMAX 

11 FORMATC ENTER NUMBER OF TEXTURES TO PROCESS .LIMIT OF',13, 

READ(*,3) T 

IF((T.LT.2) .OR. (T.GT.TMAX)) GO TO 10 
C 

12 WRITE(*,13) IK 

13 FORMATC enter THE NUMBER OF IMAGE CHANNELS .LIMIT OF',' ',I3,$) 

READ(*,3) K ' • ' 

IF((K.LT.l) .OR. (K.GT.IK)) GO TO 12 

GET ALL CHANNELS OF THE IMAGE 

DO 19 I = 1 , K 
WRITE (*,20) I 

20 FORMATC ENTER FILENAME OF THE IMAGE CHANNEL NUMBER',' ',I2,$) 

READ (*,21) IMAGE (I) 

21 FORMAT(ASO) 

OPEN(UNIT=l,FILE=IMAGE(I) ,STATUS='OLD' ,ACCESS= ' DIRECT ' ) 

DO 23 ROW = 1 . N 

READ(l'ROW) (BINPUT(ROW,COL,I) ,COL = 1 , M ) 

23 CONTINUE 

CLOSE ( UNIT = 1 ) 

19 CONTINUE 

GET THE FILTER COEFFICIENT, MEAN, AND COVARIANCE MATRICES 

DO 25 1=1, T 
WRITE(*,26) I 

26 FORMATCENTER FILENAME OF FILTER PARAMETERS SET NUMBER',' ',I2,$) 

READ(*,21) FNAME(I) 

OPEN (UNIT=2 , FILE=FNAME ( I ) , STATUS= ' OLD ' , FORM= ' FORMATTED ' ) 

DO 260 J = 1 , K 
READ(2,27) MEAN(J,I) 

27 FORMAT (FIO. 4) 

260 CONTINUE 

C 

DO 28 J=1,K 

READ(2,27) (KW( J , JJ , I) , JJ=1 ,K) 

28 CONTINUE 
C 

DO 29 PP=1,P 
DO 30 00=1,0 
DO 31 ROW=l,K 

READ(2,33) (A(ROW,COL,PP,QQ,I) ,COL=l,K) 

33 FORMAT(F10.4) 

31 CONTINUE 

30 CONTINUE 

29 CONTINUE 
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CLOSE (UNIT=2) 

32 FORMAT(3(F10.4,3X)) 

C 

25 CONTINUE 

GET THE OUTPUT FILENAMES OF THE ML AND MAP SEGMENTATION RESULTS. 

READ(*,220) MLTEST 
READ (*,220) MAPTEST 
220 FORMAT (A80) 

CONVERT A 2-D BYTE INPUT IMAGE DATA ARRAY IN THE RANGE OF -128 TO 
127 THAT REPRESENTS APPROPRIATE INTENSITY LEVELS IN THE RANGE OF 
0 TO 255 . 

DO 35 J = 1 , K 
DO 36 ROW = 1 , N 

DO 37 COL = 1 , M 

TEMP = B INPUT (ROW, COL, J) 

IF (TEMP. LT. 0.0) THEN 

TEMP = TEMP+256 

END IF 



TEXTURI 


[ROW,COL,J,l) 


= TEMP 


- MEAN! 


[J,l) 


TEXTURI 


[ROW, COL, J, 2) 


= TEMP 


- MEANI 


U.2) 



37 CONTINUE 

36 CONTINUE 

35 CONTINUE 

CALCULATION OF ERROR ESTIMATE FOR TWO TEXTURES 



46 

45 

44 

43 

421 

42 

41 

40 

C 



48 

47 

C 



C 



DO 40 I = 1 , T 

DO 41 L = 1 , N 
DO 42 LL = 1 , M 
DO 421 KK = 1 , K 

ERROR(L,LL,KK,I) = 0.0 
DO 43 III = 1 , P 
J = 1 - III 
DO 44 LLL = 1 , Q 
JJJ = LL - LLL 
DO 45 II = 1 , K 
DO 46 JJ = 1 , K 
IF(J.LE.O) J=1 
if(jjj.le.o) JJJ=1 

ERROR(L,LL,KK,I)=ERROR(L,LL,KK,I)+A(II, JJ,III,LLL,I)*TEXTUR(J, JJJ,KK,I) 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

DO 47 JJ=1,K 
DO 48 LL=1,K 



KW1(JJ,LL) 


= KW< 


;jJ,LL,l) 


KW2(jJ,LL) 


= KW( 


[JJ,LL,2) 



CONTINUE 

CONTINUE 

D1 = 1.0 

CALL LINV3F(KW1,6,1,K,K,D1,D2,WKAREA,IER) 
DETl = D1 * 2**D2 
LNU) = ALOG(DETl) 

D1 = 1.0 

CALL LINV3F(KW2,6,1,K,K,D1,DD2,WKAREA,IER) 
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DET2 = D1 * 2**DD2 
LN(2) = AL0G(DET2) 

CALCULATION OF MAXIMUM-LIKELIHOOD IMAGE SEGMENTATION 

OPEN (UNIT=3,FILE=MLTEST,STATUS=' NEW ,ACCESS= ' DIRECT ' , 

* RECL=(IM/4hMAKREC=IN) 

DO 50 ROW =1,N 
DO 51 COL = 1,M 
DO 52 I = 1,K 

EKW1(1,I) =0.0 
EKW2(1,I) = 0.0 
DO 53 J=1,K 

EKWl ( 1 , I )=EKW1 ( 1 , I )+ERROR(ROW , COL , J , 1 ) *KW1 ( J , I ) 
EKW2(1,I)=EKW2(1,I)+ERR0R(R0W,C0L,J,2)*KW2(J,I) 

53 CONTINUE 

52 CONTINUE 

PMLll =0.0 
PML22 =0.0 
DO 54 L=1,K 

PML 1 1=PML 1 1+EKWl ( 1 , L ) *ERR0R ( ROW , COL , L , 1 ) 
PML22=PML22+EKW2 ( 1 , L ) *ERR0R ( ROW , COL , L , 2 ) 

54 CONTINUE 

. PMLl =0.0 
PML2 =0.0 
PMLl = PMLll + LN(1) 

PML2 = PML22 + LN(2) 

ML(ROW,COL)=0 

PML (ROW, COL) = PML2 - PMLl 

IF (PMLl .GT. PML2) THEN 
ML(ROW,COL)= -1 
END IF 

MLI(ROW,COL) = ML(ROW,COL) 

51 CONTINUE 

WRITE(3'ROW) (ML(ROW,COL) ,COL=l,M) 

50 CONTINUE 

CLOSE (UNIT=3) 

MAXIMUM A POSTERIORI IMAGE SEGMENTATION 

0PEN(UNIT=4,FILE=MAPTEST,STATUS='NEW ,ACCESS= ' DIRECT ' , 

* RECL= (IM/4), MAXREC = IN) 

KS =10.0 
COl = 0 
CIO = 0 

DO 600 II = 1 , 5 
DO 60 1=1, N 

DO 61 J = 1 , M 
SUMl =0.0 
SUM2 =0.0 
LI = I - 3 
HI = I + 3 
LJ = J - 3 
HJ = J + 3 

IF(LI.EQ.O) LI = 1 
IF(HI.GT.N) HI = N 
IF(LJ.EQ.O) LJ = 1 
IF(HJ.GT.M) HJ = M 

AREA = (HI - LI + 1) * (HJ - LJ + 1) 

DO 62 ROW = LI , HI 
DO 63 COL = LJ , HJ 
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SUMl = SUMl - MLI(ROW,COL) 



63 CONTINUE 

62 CONTINUE 

SUMl = SUMl + MLI(I,J) 

MAP(I,J)= 0 
C 

SUM2 = PML(I,J) - ((KS/AREA)*(2*SUMl-area+l)) 
IF(SUM2,LT,0.0) THEN 
MAP(I,J) = -1 
END IF 

MLI(I,J) = MAP(I,J) 

C 

61 CONTINUE 

60 CONTINUE 
600 CONTINUE 
C 

DO 70 I = 1 , N 
DO 71 J = 1 , M 

IF(ML(I,J) .EQ. 0) THEN 

IF(MAP(I,J) .NE. ML(I,J)) COl = COl + 1 
ELSE 

IF(MAP(I,J) .NE. ML(I,J)) C10=C10+1 
END IF 

71 CONTINUE 

WRITE(4'I) (MAP(I,J),J=1,M) 

■ 70 CONTINUE 

CLOSE (UNIT=4) 

WRITER*, 64) COl, CIO 

64 FORMATC COl: M5,5 x, 'CIO: ' ,15) 

STOP 

END 
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APPENDIX E 

PROGRAMS FOR K-L TRANSFORMATION, AND ONE-CHANNEL SEGMENTATIO> 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 



c 

c 

c 



100 

101 

c 



102 



c 

c 

c 

c 



2 



c 

c 

c 

c 
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******************************************************************** 



* * 

* PROGRAM KRLV * 

* * 

* PURPOSE To implement Karhunen-Loeve transformation from a 2-D * 

* color image which size is 128 by 128 pixels. * 

* REQUIRED IMSL ROUTINES * 

* EIGRS, EQRT2S, EHOBKS , EHOUSS , UERTST, USPKD, UGETIO * 

* IMPLEMENTED BY LTJG TIMUR KUPELI Dec 1986 * 

* * 



******************************************************************** 



INTEGER I,J,L,K,N1,N2,NN,R0W,C0L,ITEMP,QC(128,128), 
MAXVAL , MINVAL , ID IF , INTVAL , N 



REAL R(l:3,l:3) ,E (1 :3 , 1 :3) , FINPUT( 1 :128 , 1 : 128 , 1 :3) , 

D(3) ,WK(3) ,TEMP,Q(1;128, 1:128, 1:3) , SLOPE 

CHARACTER*50 IMAGE (1:3) , •FNAME(1:3) 

BYTE BINPUT(1:128,1:128),QQ(1: 128, 1:128, 1:3) 

TYPE 100 

FORMAT ( 'ENTER THE NUMBER OF ROWS, COLUMNS IN THE IMAGE ',$) 

READ 101, N 
FORMAT (13) 

TYPE 102 

FORMAT (' ENTER THE NUMBER OF CHANNELS IN THE IMAGE ',$) 

READ 101, K • 

GET THE FILENAMES OF THE RED , GREEN, AND BLUE COMPONENTS 
OF THE COLOR IMAGE. 

DO 1 I = 1 , K 
WRITE(*,2) I 

FORMAT (' ENTER THE FILENAME OF THE IMAGE CHANNEL NUMBER ',I2,$) 
READ(^,3) IMAGE(I) 

WRITE(*,3) IMAGE(I) 

WRITE (*,45) 

FORMAT (A50) 

CONVERT THE IMAGE BTYE FORMAT TO THE REAL NUMBER FORMAT 
OPEN(UNIT=l, FILE=IMAGE(I) , STATUS= ' OLD ' , ACCESS= ' DIRECT ' ) 

DO 5 ROW = 1 , N 

READ(l'ROW) (BINPUT(ROW,COL) ,COL=l,N) 

DO 6 COL = 1 , N 

TEMP = B INPUT (ROW, COL) 

IF (TEMP. LT. 0.0) THEN 

TEMP = TEMP +256 

END IF 

F INPUT (ROW, COL, I) = TEMP 
CONTINUE 

CONTINUE 
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c 



23 

22 

21 

20 



39 

30 

31 



CLOSE (UNIT = 1) 

CONTINUE 

CALCULATE THE CORRELATION MATRIX 

NN = N * N 
DO 20 I = 1 , K 

DO 21 J = 1 , K 
R(I,J) = 0 
DO 22 N1 = 1 , N 

DO 23 N2 = 1 , N 

R(I,J) = R(I,J) + FINPUT(N1,N2,I)*FINPUT(N1,N2,J) 
CONTINUE 
CONTINUE 

R(I,J) = R(I,J) / NN 
CONTINUE 
CONTINUE 

WRITE (*,45) 

WRITE(*,30) 

WRITE (*,39) 

FORMATC ') 

FORMAT (‘ ',5X,' THE CORRELATION MATRIX' , $) 

WRITE/*, 3i; ((R(I,J),J=1,K),I=1,K) 

FORMAT ( <K>(F9.2,4X)) 

WRITE (*,39) 

CALCULATE EIGENVALUES AND EIGENVECTORS OF THE 
CORRELATION MATRIX 

JOHN = 11 

CALL EIGRS(R,K, JOBN,D,E,K,WK,IER) 

SORT THE EIGENVALUES IN DECREASING ORDER 

TEMP = D(l) 

D(l) = D(3) 

D(3) = TEMP • 

DO 49 I = 1 , K 
TEMP 



49 


CONTINUE 




WRITE(*,39) 

WRITE(*,45) 


45 


FORMATC ') 




WRITE(*,41) 
WRITE (*,39J 


41 


FORMATC ‘,5X 




DO 46 J=1,K 
WRITE (*,42) 


42 


FORMAT(5X,F9. 


46 


CONTINUE 



1 , 1 ) = 
1,3) = 



= E 1,1 
= E(I,3) 
TEMP 



' THE EIGENVALUES 
P(J) 



,$) 



WRITE(*,39) 

WRITE (*,45) 

WRITE (*,43) 

WRITE(*,39j 

43 FORMATC ',5X,' THE EIGENVECTORS', $) 

WRITE (*.44) ((E(I,J) ,J=1,K) ,I=1,K) 

44 FORMAT(3(^F9.2,4X)) 

WRITE (*,39) 



USE THE TRANSPOSE OF THE EIGENVECTORS MATRIX TO IMPLEMENT 
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KARHUNEN-LOEVE TRANSFORMATION . 

DO 50 I = 1 , K 

DO 51 N1 = 1 , N 

DO 52 N2 = 1 , N 

Q(N1,N2,I) = 0.0 
L = 4 - I 
DO 53 J = 1 , K 

Q(N1,N2,I) = Q(N1,N2,I) + E ( J , I )*FINPUT(N1 ,N2 , J) 
53 CONTINUE 

52 CONTINUE 

51 CONTINUE 

50 CONTINUE 

DO 60 I = 1 , K 
WRITE(*,61) I 

61 FORMATC ENTER FILENAME OF THE TRANSFORMED IMAGE M2,$) 

READ(*,3) FNAME(I) 

WRITE(*,3) FNAME(I) 

WRITE(*,45) 

CONVERT THE TRANSFORMED IMAGE TO BYTE FORMAT 



DO 62 N1 = 1 , N 
DO 63 N2 = 1 , N 

QC(N1,N2) = JNINT(Q(N1,N2,D) 

63 CONTINUE 

62 CONTINUE 

SCALE THE TRANSFORMED IMAGE TO BE WITHIN DISPLAY RANGE 



IF. (I .EQ. 1) THEN 
MAXVAL = QC(1,1) 

MINVAL = QC(1,1) 

DO ROW = 1 , N 
DO COL = 1 , N 

IF (QC(ROW,COL) .GT. MAXVAL) THEN 
MAXVAL = QC(ROW,COL) 

ELSE IF (QC(ROW,COL) .LT. MINVAL) THEN 
MINVAL = QC (ROW, COL) 

END IF 
ENDDO 
ENDDO 

INTVAL = MAXVAL - MINVAL 
SLOPE = 255 / REAL (INTVAL) 

DO ROW = 1 , N 
DO COL = 1 , N 

IF (QC(ROW,COL) .EQ. MINVAL ) THEN 
QC(ROW,C(^) = 0 

ELSE IF (QC(ROW,COL) .EQ. MAXVAL) THEN 
QC(ROW,COL) = 255 

ELSE 



IDIF = QC(ROW,COL) - MINVAL 
QC(ROW,COL) = INT(SLOPE*IDIF) 

END IF 
ENDDO 
ENDDO 
ELSE 

MINVAL = QC(1,1) 

DO ROW = 1 , N 
DO COL = 1 , N 

IF (QC(ROW,COL) .LT. MINVAL) THEN 
MINVAL = QC (ROW, COL) 

END IF 
ENDDO 
ENDDO 

DO ROW = 1 , N 
DO COL = 1 , N 

IDIF = QC (ROW, COL) - MINVAL 
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QC(ROW,COL) = INT(SLOPE*IDIF) 



ENDDO 
ENDDO 
END IF 
C 

DO 64 N1 = 1 , N 
DO 65 N2 = 1 , N 

ITEMP = QC(N1,N2) 

IF (ITEMP .GT.127) THEN 

ITEMP = ITEMP - 256 

END IF 

QQ(N1,N2,I) = ITEMP 

65 CONTINUE 

64 CONTINUE 

C 

0PEN(UNIT=2,FILE=FNAME(I) ,STATUS='NEW ,ACCESS= ' DIRECT ' , 
" RECL= (N/4), MAXREC=N ) 



DO 66 N1 = 1 , N 
WRITE (2 ‘Nl) (QQ(N1,N2,I),N2=1,N) 
66 CONTINUE 

C 

CLOSE (UNIT= 2) 

C 

60 CONTINUE 

C 

STOP 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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c 
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PROGRAM FLTR2 

PURPOSE To develop two sets of filter parameters from a 2-D single 
channel image which size is 128 by 128 pixels. These para- 
meters are the mean, the covariance, and the filter 
coefficients . 

REQUIRED IMSL ROUTINES 

LEQT2F, LUDATN, LUELMN, LUREFN 

IMPLEMENTED BY LTJG TIMUR KUPELI Sep 1986 



* 
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INTEGER*4 P,Q,R0W,C0L,STARTN(2) ,STARTM(2) ,ENDN(2) ,ENDM(2) , IDGT, 

PQ , RNROW , RNCOL , RROW , CCOL , X , Y , 

HNNO , HMMO , HCOL , HROW , LNNO , LNMO , LCOL , LROW , RN , RM , lER 



* 

k 



REAL*4 SUM , TEMP , F INPUT (128,128), MEAN ( 2 ) , WKAREA ( 28 ) , 

R(4,4),SI(4,1) ,IMATRX(1,1) ,B00, 
KW(1,1),A(4,1),ISIZE 



CHARACTER*50 IMAGE,FILTER(2) 



BYTE BINPUT(128) 



T = 2 



GET PROGRAM INPUT PARAMETERS. 

TYPE 10 

10 FORMATC enter number of filter rows from 2-4 DESIRED :',$) 

READ 11, P 

11 FORMAT(I3) 

TYPE 12 

12 FORMATC enter NUMBER OF FILTER COLUMNS FROM 2-4 DESIRED :',$) 

READ 11, Q 

TYPE 13 

13 FORMATC ENTER NUMBER OF ROWS IN IMAGE ',$) 

READ 11, N 

TYPE 14 

14 FORMATC ENTER NUMBER OF COLUMNS IN IMAGE ',$) 

READ 11, M 

GET FILENAME OF SINGLE CHANNEL IMAGE 
TYPE 15 

15 FORMATC ENTER FILENAME OF IMAGE ',$) 

READ 16, IMAGE 

16 FORMAT(A50) 

CONVERT THE IMAGE BTYE FORMAT TO REAL NUMBER FORMAT 

OPEN(UNIT=l,FILE=IMAGE,STATUS='OLD' , ACCESS = 'DIRECT') 

DO 17 ROW = 1 , N 

READ (I'ROW) (BINPUT(COL),COL=l,M) 

DO 18 COL = 1 , M 

TEMP = BINPUT(COL) 

IF (TEMP. LT. 0.0) THEN 
TEMP = TEMP + 256 
END IF 

FINPUT(ROW,COL) = TEMP 
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18 

17 

20 

21 

22 

23 

24 

26 

25 

27 

29 

28 

291 



CONTINUE 

CONTINUE 
CLOSE (UNIT=1) 

GET 'T' AREAS FOR WHICH FILTERS ARE DESIRED AND OUTPUT FILENAME FOR 
EACH AREA'S FILTER COEFFICIENTS AND COVARIANCE MATRIX. 

DO 19 I = 1,T 
WRITE (*,20) I 

FORMATC ENTER UPPER-LEFT ROW FOR AREA ',12,':',$) 

READ(*,11) STARTN(I) 

WRITE (*,21) I 

FORMAT ( ' ENTER UPPER-LEFT COLUMN FOR AREA' ,12, ' : ' ,$) 

READ (*,11) STARTM(I) 

WRITE (*,22) I 

FORMAT ( ' ENTER LOWER-RIGHT ROW FOR AREA' ,12, ' ; ' ,$) 

READ(*,11) ENDN(I) 

WRITE (*,23) I 

FORMAT ( ' ENTER LOWER-RIGTH COLUMN FOR AREA' ,12, ' : ' ,$) 

READ(*,11) ENDM(I) 

WRITE (*,24) I 

FORMATC ENTER OUTPUT FILE-SPEC FOR FILTER' ,12, ';' ,$) 

READ(*,16) FILTER(I) 

FIND THE MEAN VECTOR OF ESTIMATION WINDOW AREA OF IMAGE 

ISIZE = (ENDN(I) - STARTN(I) + 1)*(ENDM(I) - STARTM(I) + 1) 

SUM = 0.0 

DO 25 L = STARTN(I) ,ENDN(I) 

DO 26 J = STARTM(I) ,ENDM(I) 

SUM = SUM + FINPUT(L,J) 

CONTINUE 

CONTINUE 

MEAN (I) = SUM / ISIZE 
WRITE(*,27) I,MEAN(I) 

FORMAT(-' MEAN' , ' ( ' ,12, ')',':' ,F9.2,$) 

CORRECT THE IMAGE TO BE ZERO MEAN 

DO 28 L = STARTN(I), ENDN(I) 

DO 29 J = STARTM(I) ,ENDM(I) 

FINPUT(L,J) = FINPUT(L,J) - MEAN(I) 

CONTINUE 

CONTINUE 

DETERMINE CORRELATION MATRIX 
WRITE (*,291) I 

FORMATC' CORRELATION. MATRIX' ,12, $) 

PQ = P * Q 
DO 30 RNROW = 1,P 
X = Q * (RNROW-1) 

DO 31 RNCOL = 1,P 

NO = RNROW - RNCOL 
Y = Q * (RNCOL-1) 

DO 32 RMROW = 1,Q 
RN = X + RMROW 



3 RMCOL =1,0 
MO = RMROW - RMCOL 


RM = 
LNNO 


Y + RMCOL 
= STARTN(] 


[) 


+ NO 


LMMO 


= STARTM(] 


[) 


+ MO 


HNNO 


= ENDN(I) 


+ 


NO 


HMMO 


= ENDM(I) 


+ 


MO 
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35 

34 



33 

32 

31 

30 



LCOL =MAX0(STARTM(I),LMM0) 

HCOL =MIN0(ENDM(I) ,HMM0) 

LROW =MAX0(STARTN(I) ,LNN0) 

HROW =MIN0(ENDN(I),HNN0) 

SUM = 0.0 

DO 34 ROW = LROW, HROW 

NNO = ROW - NO 
DO 35 COL = LCOL, HCOL 
MMO = COL - MO 

SUM = SUM + FINPUT(ROW,COL)*FINPUT(NNO,MMO) 

CONTINUE 

CONTINUE 

R(RN,RM) = SUM / ISIZE 



CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 



RESET SI(J,1) TO HAVE 1 IN THE FIRST ROW AND 0 IN ALL OTHERS 

DO 41 J = 1,PQ 

IF (J.EQ.l) THEN 

SI(J,1) = 1.0 

IMATRXU.l) = 1-0 

ELSE 

SI(J,1) = 0.0 

END IF 

41 CONTINUE 



THE FORMATS BELOW MUST BE MODIFIED TO USE DIFFERENT FILTER 
SIZE THAN 2 by 2 PIXELS.- 

IDGT = 1 
DO 36 K=1,PQ 

WRITE(*,37) (R(K,J),J=1,PQ) 

WRITE(*,38) 

37 FORMATC ',4F9.2,$) 

38 FORMAT(' ',$) 

36 CONTINUE 



360 



SOLVE EQUATION [ R ] * [ B ] = [ SI ] 

. CALL LEQT2F (R, IDGT , PQ , PQ ,SI , 3 , WKAREA, lER) 

WRITE(*,360) (SI(J,1) , J=1,PQ) 

FORMAT (FI 0.2) 

SOLVE EQUATION BOO * KW = I 
BOO = SI(1,1) 

GET THE FILTER COEFFICIENTS USING THE EQUATION 
[ A ]= [ B ] * KW 

KW(1,1) = IMATRX(1,1) / BOO 
WRITE (^,*) 'COVARIANCE' 

WRITE(*,360) KW(1,1) 



DO 43 ROW = 1,PQ 

A(ROW,l) = SI(R0W,1) * KW(1,1) 
43 CONTINUE 



WRITE OUT MEAN, COVARIANCE, AND FILTER COEFFICIENTS TO THE 
USER INPUT FILE. 

OPEN (UNIT=2,FILE=FILTER(I) ,STATUS='NEW' ,CARRIAGECONTROL='LIST' , 
FORM= ' FORMATTED ' ) 
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500 



51 

52 
C 

C 



53 

54 
19 



WRITE(2,50) KW(1,1) 
FORMAT ( F10.5) 
WRITE(2,500j MEAN(I) 
FORMAT ( F10.5) 

DO 51 J = 1,PQ 

WRITE(2,52) A(J,1) 
CONTINUE 
FORMAT (F10.5) 



CLOSE (UNIT=2) ' 

WRITE (*,53) I 

FORMAT FILTER COEFFICIENTS :',I2,$) 
WRITE(*,54) (A(J,1),J=1,PQ) 

FORMAT (FIO. 5) 

CONTINUE 

STOP 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 



c 

c 

c 
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c 
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C 

6 
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C 
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C 

C 

c 
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* PROGRAM SGMT2 

* 

* PURPOSE To segment a single image with two textures given the 

* image, the image dimensions, the filter dimensions, 

* and two sets or filter parameters to be used. 

* REQUIRED IMSL ROUTINES 

* NONE 

* 

* IMPLEMENTED BY LTJG TIMUR KUPELI Sep 1986 

* 
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INTEGER IN,N,IM,Pl,P,Ql,Q,TMAX,T,PQ,ROW,COL,I,J,JJ, JJJ,L,LL, 

LLL , LLLL , COUNT , LI , HI , L J , HJ 

REAL KW(2) ,MEAN(2),AA(4,1,2) ,TEMP, SUMl ,SUM2 ,TEXTUR(128 , 128, 2) , 
LN(2) ,ERROR(128, 128,2), PML1,PML2,PML(128, 128) , 
AREA,KS,AU/2) 

CHARACTER*50 IMAGE , FNAME ,MLTEST ,MAPTEST . 

BYTE BINPUT(128,128) ,ML(128, 128) ,MAP(128 , 128) ,MLI(128, 128) 

IN = 128 
IM = 128 
PI = 4 
Q1 = 4 
TMAX = 2 



GET THE INPUT PARAMETERS OF THE PROGRAM 
WRITE(*,2) IN 

FORMATC ENTER THE NUMBER OF ROWS IN IMAGE. LIMIT OF ' , 13 , ' ; ' , $) 
READ(*,3) N 
FORMAT (13) 

IF((N.LT.l) .OR. (N.GT.IN)) GOTO 1 
WRITE(*,5) IM 

FORMATC enter THE NUMBER OF COLUMNS IN IMAGE. LIMIT OF ', 13 ,';',$ ) 
READ(*,3) M 

IF((M.LT.l) .OR. (M.GT.IM)) GOTO 4 
WRITE(*,7) PI 

FORMATC ENTER THE NUMBER OF ROWS IN FILTER. LIMIT OF ' , 13 , ' ; ' , $) 
READ(*,3) P 

IF((P.LT.2) .OR. (P.GT.Pl)) GOTO 6 
WRITE(*,9) Q1 

FORMATC ENTER THE NUMBER OF COLUMNS IN FILTER. LIMIT OF ' , 13 , ' : ' , $) 
READ(*,3) Q 

IF((Q.LT.2) .OR. (Q.GT.Ql)) GOTO 8 
WRITE (*,11) TMAX 

FORMATC ENTER NUMBER OF TEXTURES TO PROCESS .LIMIT OF ' , 13 , ' : ' , $) 
READ(*,3) T 

IF((T.LT.2) .OR. (T.GT.TMAX)) GO TO 10 
GET THE SINGLE -CHANNEL IMAGE 



PQ = p * 

WRITE (*,26) 

FORMATC enter FILENAME OF IMAGE 
READ (*,21) IMAGE 
FORMAT (A50) 



,$) 
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c 



c 

22 

C 



0PEN(UNIT=1,FILE=IMAGE,STATUS='0LD' ,ACCESS= ‘ DIRECT ' ) 
DO 22 ROW = 1 . N 

READ(l'ROW) (BINPUT(ROW,COL),COL = 1 , M ) 

CONTINUE 

CLOSE ( UNIT = 1 ) 

GET THE FILTER COEFFICIENTS , MEANS , AND COVARIANCES 



DO 23 I = 1 , T 

WRITE (*,24) I 

24 FORMAT (' ENTER THE FILENAME OF FILTER PARAMETERS SET NUMBER',12, ' : ' 

READ (*,21) FNAME 

0PEN(UNIT=2 , FILE=FNAME , STATUS= ' OLD ' ) 

READ(2,240) KW(I) 

240 FORMAT(FIO.S) 

READ (2, 240) MEAN(I) 

DO 242 KK= 1,PQ 
READ(2,240) AA(KK,1,I) 

242 CONTINUE 

CLOSE (UNIT = 2) 

WRITE (*,241) KW(I) ,MEAN(I),AA(1,1,I) , AA(2 , 1 , I ) , AA(3 , 1 , I ) , AA(4, 1 , I ) 

241 FORMAT(6F10.5) 

23 CONTINUE 

READ (*,220) MLTEST 
READ (*,220) MAPTEST 
220 FORMAT (A80) 

CONVERT A 2-D BYTE INPUT IMAGE DATA IN THE RANGE OF -128 TO 127 THAT 
REPRESENTS APPROPRIATE INTENSITY LEVELS IN THE RANGE OF 0 TO 255 



DO 30 ROW = 1 , N 

DO 31 COL = 1 , M 

TEMP = B INPUT (ROW, COL) 

IF (TEMP. LT. 0.0) TEMP = TEMP+256 
TEXTUR(ROW,COL,l) = TEMP - MEAN(l) 
TEXTUR(ROW,COL,2) = TEMP - MEAN(2) 

31 CONTINUE 

30 CONTINUE 



CALCULATOIN OF ERROR ESTIMATE FOR TWO TEXTURES 

DO 40 I = 1 , T 
J = 1 

DO 41 JJ = 1 , P 

DO 42 JJJ = 1 , Q 

A(JJ,JJJ) = AA(J,1,I) 

J = J + 1 
42 CONTINUE 

41 CONTINUE 

C 

DO 43 L = 1 , N 
DO 44 LL = 1 , M 

ERROR(L,LL,I) = 0.0 
DO 45 LLL = 1 , P 
J = L - LLL 
DO 46 LLLL =1,0 

JJ = LL - LLLL 

IF((J.GT.O) .AND. (JJ.GT.O)) THEN 
ERROR(L,LL,I) =ERROR(L,LL, I)+ A(LLL,LLLL)*TEXTUR( J, JJ, I) 
END IF 

46 CONTINUE 
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45 CONTINUE 

44 CONTINUE 

43 CONTINUE 

C 

LN(I) = ALOG(KW(D) 
40 CONTINUE 



CALCULATION OF MAXIMUM LIKELIHOOD IMAGE SEGMENTATION 



51 

50 



0PEN(UNIT=3,FILE=MLTEST,STATUS='NEW ,ACCESS=' DIRECT' , 
RECL=(IM/4) ,MAXREC=IN) 

L = 1 
LL = 2 

DO 50 I = 1 , N 
DO 51 J = 1 , M 
PMLl =0.0 
PML2 =0.0 

PMLl = ((ERR0R(I,J,L)**2) 

PML2 = ((ERR0R(I, J,LL)**2 
PML(I,J) = PML2 - PMLl 
IF (PMLl .GT. PML2) THEN 
ML(I,J)= -1 
END IF 

MLI(I,J) = ML(I,J) 

CONTINUE 

WRITE(3'I) (ML(I,J) ,J=1,M) 

CONTINUE 



'KW(LL)) +^Ln1lL) 



CLOSE (UNIT=3) 



MAXIMUM A POSTERIORI IMAGE SEGMENTATION 



C 



C 



C 

C 
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C 



OPEN (UNIT=4,FILE=MAPTEST,STATUS=' NEW ,ACCESS= ' DIRECT ' , 
RECL= (IM/4), MAXREC = IN) 

KS = 100 

DO 600 II = 1 , 5 

DO 60 I = 1 ,N 

DO 61 J = 1 , M 
SUMl =0.0 
SUM2 =0.0 
LI = I - 3 

HI = I + 3 

LJ = J - 3 

HJ = J + 3 



IF(LI.EQ.O) LI = 1 
IF(HI.GT.N) HI = N 
IF(LJ.EQ.O) LJ = 1 
IF(HJ.GT.M) HJ = M 

AREA = (HI - LI + 1) * (HJ - LJ +1) 

DO 62 ROW = LI , HI 
DO 63 COL = LJ , HJ 

IF((R0W.NE.I) .OR. (COL.NE.J)) THEN 
SUMl = SUMl - MLI(R0W,C0L) 

END IF 
CONTINUE 
CONTINUE 

MAP(I,J)=0 

SUM2 = PML(I.J) - ((KS/AREA)*(2*SUM1-AREA+D) 
IF(SUM2.LT.O) THEN 
MAP(I,J) = -1 
END IF 
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MLI(I,J) = MAP(I,J) 

61 CONTINUE 

60 CONTINUE 

600 CONTINUE 

C 

DO 70 I = 1 , N 

WRITE(4'I) (MAP(I,J) ,J=1,M) 

70 CONTINUE 
STOP 
END 
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